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Accurate models of alkali and halide ions in aqueous solution are necessary for computer simulations of a 
broad variety of systems. Previous efforts to develop ion force fields have generally focused on reproducing 
experimental measurements of aqueous solution properties such as hydration free energies and ion-water 
distribution functions. This dependency limits transferability of the resulting parameters because of the 
variety and known limitations of water models. We present a solvent-independent approach to calibrating ion 
parameters based exclusively on crystal lattice properties. Our procedure relies on minimization of lattice 
sums to calculate lattice energies and interionic distances instead of equilibrium ensemble simulations of dense 
fluids. The gain in computational efficiency enables simultaneous optimization of all parameters for Li+, 
Na+, K+, Rb+, Cs+, F-, C1-, Br-, and I- subject to constraints that enforce consistency with periodic table 
trends. We demonstrate the method by presenting lattice-derived parameters for the primitive model and the 
Lennard- Jones model with Lorentz-Bcrthelot mixing rules. The resulting parameters successfully reproduce 
the lattice properties used to derive them and are free from the influence of any water model. To assess 
the transferability of the Lennard- Jones parameters to aqueous systems, we used them to estimate hydration 
free energies and found that the results were in quantitative agreement with experimentally measured values. 
These lattice-derived parameters are applicable in simulations where coupling of ion parameters to a particular 
solvent model is undesirable. The simplicity and low computational demands of the calibration procedure 
make it suitable for parametrization of crystallizable ions in a variety of force fields. 
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I. INTRODUCTION 

Alkali and halide ions play important roles in biologi- 
cal and physico-chemical systems that include proteirfQHl 
nucleic acicPKin lipicfCLJ an( j carbohydrate^ solutionis, 

salt crystal^!, molten salts^, electrolyte^, and liquid- 
vapor interfaces^. Computer simulations are useful for 
developing a molecular scale description and understand- 
ing of electrolyte dependencies and ion-mediated inter- 
actions in these systems. Most classical simulation ap- 
proaches to modeling alkali and halide ions employ the 
Born-Oppenheimer approximation where the ions are 
hard spheres or van der Waals spheres with a charge of 
±e. van der Waals interactions are commonly modeled 
using the empirical Lennard-Jones 12-6 potential. These 



a 'Electronic mail: albert.mao@gmail.com 

b ) Electronic mail: fpappu@wustl.edu| PEone: (314) 935-7958 



models for ions are used with either explicit!^ or contin- 
uurr ji9H22] (i m pii c it) descriptions of the surrounding sol- 
vent. In order to achieve accuracy in simulation results, 
one needs reliable hard sphere or van der Waals param- 
eters. 

Experiments that measure structural and excess ther- 
modynamic properties of electrolyte solutions can pro- 
vide constraints for the calibration of these parameters. 
Numerous collections of parameters have been developed 
for ions in aqueous solutions where the relevant con- 
straints come from gas phase ion-water binding energies 
and geometries^2tt26|, hydration free energies^SU anc j en _ 
tropieiP^, structur al prop erties such as water-ion pair dis- 
tribution function d 25 l 27 tt 221 anc l transport properties re- 
garding the degree of hydration and ion association^. 
The resulting parameters are intricately dependent on 
the water model used in the parameterization procedure. 
This dependence is awkward because water models them- 
selves are more complicated and require more parameters 
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than any one alkali or halide ion. The need for transfer- 
able and generally applicable parameters is prominent in 
applications such as biomolccular simulation, where mat- 
ter besides water and ions is present. Given the sheer 
number, diversity, and known limitations^ of available 
water models, it is difficult to be confident that ion pa- 
rameters derived using a particular water model reflect 
the intrinsic properties of the ions that are also transfer- 
able for use in a specific simulation system. 

The interionic distance and lattice energy of alkali 
halide salt crystals are properties that do not require 
any consideration of the specific model used for solvent 
molecules. They constitute a set of measurable observ- 
ables that constrain the length and energy scales for 
van der Waals interactions of non-polarizable alkali and 
halide ions. In this work, we simultaneously obtain values 
for the sphere diameter (a) and well depth (e) parame- 
ters of five alkali cations (Li + , Na + , K + , Rb + , Cs + ) and 
four halide anions (F - , CP, Br~, I - ). These param- 
eters, which are designed for use with a hard sphere or 
Lennard- Jones 12-6 potential based on Lorentz-Berthelot 
mixing rules, were obtained using lattice properties as the 
only calibration targets and are therefore independent of 
any specific water or solvent model. Our calibration pro- 
cedure relies on minimization of lattice sums to compute 
the lattice energies and interionic distances. It requires 
modest computational resources compared to approaches 
involving explicit construction and simulation of periodic 
crystals or dense fluids. We assess the transferability of 
the derived parameters by determining minimum energy 
lattice configurations and testing the accuracy of single- 
ion hydration free energies across three water models es- 
timated using bicubic surfaces constructed by Joung and 
CheatharrP! 



II. METHODS 

The fitting is accomplished through minimization of a 
calibration objective function that maps any candidate 
parameter set to one real number quantifying deviation 
from experimental measurements of lattice observables. 
Each evaluation of the objective function itself involves 
minimization of every salt's parameter-dependent poten- 
tial energy to calculate its lattice energy and interionic 
distance. Since this minimization does not account for 
thermal fluctuations, it is a ground state calculation suit- 
able for comparison with experimental lattice energies 
and interionic distances measured at absolute zero. The 
following sections describe these steps in detail. 



A. Calibration targets 

Experimental data for the twenty alkali halide salts 
arising from combinations of Li + , Na + , K + , Rb + , or 
Cs + with F - , CP, Br~, or I form the basis of our 
calibration. All of these salts form cubic crystals whose 



structures are described by one interionic distance (ID), 
which is defined as the distance between centers of two 
nearest-neighbor ions of opposite charge. The cations 
and anions are arranged in interpenetrating simple cu- 
bic lattices (denoted BCC, since the unit cell is body- 
centered) for CsCl, CsBr, and Csl and interpenetrating 
face centered cubic (FCC) lattices for the other seventeen 
salts. Sirdeshmukh et al. have compiled X-ray diffrac- 
tion measurements of interionic distances^!, while Jenk- 
ins and Roobottom have gathered lattice energy (LE) 
measurements^ derived through the Born-Fajans-Haber 
thermochemical correlation^. Ghate^ extrapolated in- 
terionic distances to K, a reduction of ~1% relative 
to room temperature values. The lattice energies have 
even lower temperature sensitivity, changing by less than 
0.1% from room temperature to KpS. Therefore, we 
treat all parameters as temperature-independent quanti- 
ties and adopt the extrapolated interionic distances and 
room temperature lattice energies, listed in Table |TJ as 
our calibration targets because they closely approximate 
the K values that result from minimizing the potential 
energy of a crystal lattice. Since there are forty inde- 
pendent measurements and nine ions, models with four 
or fewer parameters per ion are overdetermined, which is 
a desirable characteristic in discouraging overfitting and 
promoting transferability. 

TABLE I. Lattice energy (LE) and interionic distance (ID) 
measurements used as calibration targets. For each salt, the 
top number is the negative lattice energy and the bottom 
number is the interionic distance. 



— LE (kcal/mol) 
ID (A) 


F 


CI 


Br 


I 


Li 


250.7 


206.5 


196.0 


182.6 




1.996 


2.539 


2.713 


2.951 


Na 


222.3 


188.8 


180.2 


168.5 




2.295 


2.789 


2.954 


3.194 


K 


198.1 


172.1 


165.2 


155.4 




2.648 


3.116 


3.262 


3.489 


Rb 


190.0 


166.1 


159.7 


151.1 




2.789 


3.259 


3.410 


3.628 


Cs 


181.4 


160.1 


154.6 


146.5 




2.982 


3.523 


3.668 


3.898 



B. Ion interaction models 

We focus on two common models for ions in molecular 
simulations. Both models are pairwise additive potentials 
where a short-range interaction is superimposed upon the 
Coulomb electrostatic interaction. For a pair of particles 
denoted as i and j, 

Dg«(r«) ee *p (1) 

where k = e 2 /47re ~ 332.06 kcal/mol, the valences z are 
1 for alkali cations and —1 for halide anions, and r.y is the 
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distance between their centers. The dielectric constant is 
uniformly 1 because the salt crystals are modeled in the 
absence of other matter. 

In the primitive model (PM), ions are hard spheres 
which cannot overlap. Each ion has one parameter, its 
diameter a: 



Ur™(r ij ) = U!nr ij ) + 



if 



otherwise 



(2) 



In the Lennard- Jones (LJ) model, ions exhibit short- 
range attractive van der Waals interactions that compete 
against a repulsive barrier. Each ion has two param- 
eters, a and e, which respectively describe the length 
and energy scales of its interactions. We adopt Lorentz- 
Berthelot mixing rules, which specify arithmetic means 
for a and geometric means for e, to combine the param- 
eters of two ions into one pairwise interaction. 



/ \ 12 













(3) 
(4) 



C. Lattice energy and interionic distance calculations 

Due to the symmetry of a periodic alkali halide crystal 
lattice, the contribution from one ion to the total po- 
tential energy is the same for every cation and for every 
anion. This contribution is equal to the potential energy 
of one ion in the field generated by all other ions in the 
lattice, divided by two to correct for double counting. 
The sum of the contributions from one cation and one 
anion, respectively denoted as c and a, gives the poten- 
tial energy per salt pair as a function of the tentative 
interionic distance d: 



{c,a} {lattice} 



(5) 



Proceeding in a manner analogous to the derivation of 
the classic Born-Lande equation, we compute the cali- 
bration observables by minimizing this intensive poten- 
tial energy with respect to d. The minimum energy is the 
lattice energy (LE) and optimal value of d is the interionic 
distance (ID). Since the shape of U ca (d) depends on the 
parameters "P, the calibration observables are functions 
of V . In Equations |6j [7J and 12 we make this functional 
relation explicit. V ca denotes parameters pertaining to 
the cation c and anion a and consists of the pair (cr c , a a ) 
for the primitive model and the quadruple (<j c , e c , a a , e a ) 
for the Lcnnard-Jones model. 



ID ca (V ca ) = arg mmU ca (d) 

d>0 

LEcaCPca) = min[/ ca (d) 

d>0 



(G) 



For the primitive model, since the Coulomb interac- 
tion draws the lattice together as tightly as possible, the 
optimal value of d is the one at which lattice ions come 
into contact with each other: 



oo if d < max {scr c , sa a , \{p c + cr a )} 
— ^fi- otherwise 



ID ca (cr c , a a ) = max <j sa c , sa a , -(a e + <r a ) 



1 



ID™(* Cl a a ) 



(7) 



where s is 1/^2 for FCC lattices and for BCC 

lattices and b\, the Madelung constant, is approximately 
1.7476 for FCC lattices and 1.7627 for BCC lattices. 

For the Lennard- Jones model, the intensive potential 
energy can be written as a rational function of d because 
the lattice sums bi 2 and of the Lennard- Jones potential 
are numerical constants and can be precomputed: 



hk 
~d 

{c,a} 



(8) 



Unlike the classic Madelung constant, these lattice sums 
are absolutely convergent!^. They can be approximated 
by summing over a finite cube centered at and omitting 
the origin, but care must be taken to prevent inexact 
floating point arithmetic from making the summation 
converge to an inaccurate resurP^. We use exact rational 
arithmetic for accumulating these sums and convert the 
final totals to IEEE 754 double-precision floating point 
numbers. Double-precision floating point arithmetic is 
used for all other numerical operations in this study. We 
sum over a cube with 601 3 ions, keeping separate totals 
for "even" and "odd" lattice sites, to obtain the numeri- 
cal values for FCC lattices in Equation [9j 

be,cc = b 6 . aa w 1.8067 6 6jCQ = 6 6 , ac w 6.5952 
b12.ee = b 12 . aa w 0.1896 b 12 . ca = b 12 . ac » 6.0126 (9) 

Likewise, we perform separate sums over cubes with 601 3 
and 600 3 ions for "like" and "unlike" lattice sites, respec- 
tively, to obtain the corresponding BCC lattice sums. 
Equation 
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presents these sums multiplied by \/3/2, the 
ratio of interionic distance to lattice constant, raised to 
the sixth power for be or twelfth power for b\2'. 

b 6 ,cc = b 6iaa « 3.5446 6 6 , ca = b 6iac ps 8.7091 

&12,ec = 6l2,aa « 1-1038 b l2 ,ca = &12,ac ~ 8.0103 (10) 

Decreasing cube side lengths by a factor of three changes 
the resulting sums by an amount less than 1.2 x 10 -6 . 
Minimizing (d) requires the real positive roots of 



the polynomial Q in Equation 11 which was obtained by 
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implementing the mixing rules in Equation [3] and factor- 
ing the derivative dU^ /dd: 



11 



Qca{d) = -256hkd 

-96 b 6yCay /e c e a (a c + <r a f 

+ 32 (b 6 ,ccec<J 6 c + b 6taa e a a 6 a )] d 6 

\12 



bl2, ca^/^c^a (&c H 

2048 (b 12iCC e c al 2 + b 12 , aa e a af)] 



(11) 



Given particular numerical values for the parameters cr c , 
e c , a a , and e a , lattice sums 6, and electrostatic constant 
k, Q becomes an eleventh-degree polynomial in d whose 
roots can be obtained via standard methods such as the 
Jenkins- Traub algorithm^!. In cases where the derivative 
has multiple positive roots, the smallest one is taken to 
define the interionic distance: 

ID^(a c ,e c ,a a ,e a ) = min{d € R >0 : Q ca (d) = 0} 
LE L J (a c , e c> a a , e a ) = (lD™(a c , e c , a a ,e a )) (12) 

If no positive roots exist, the situation corresponds to an 
unstable crystal with undefined calibration observables. 



D. Calibration objective function and parameter 
constraints 

The root mean square relative deviation from 
O measured , the measured values of calibration targets 
given in Table |TJ is used as the objective function whose 
minimization yields the optimized parameters. 



F(P)s 



\ 



{salts} {ID,LE} 



40 



E 



o 



C a \ Pca) 
measured 



(13) 



In Equation 13 the outer sum is over all twenty cation- 



anion salt pairs, and the inner sum is over the two ob- 
servables ID and LE. We used relative instead of abso- 
lute deviations because they put all the observables on 
an equal footing; biases due to differing magnitudes and 
units are naturally eliminated. As a special case, pa- 
rameters that cause any salt crystal to be unstable are 
defined to have an infinite calibration objective function 
value. The domain {V} of the function consists of the 
nine-dimensional space {cr} 9 for the primitive model and 
the eighteen-dimensional space {a, e} 9 for the Lennard- 
Jones model. 

Following the reasoning of Peng et alP^, we constrain 
the domain to maintain consistency with periodic table 
trends. Ions increase in size going down their respective 
groups of the periodic table, and cations are smaller than 
their isoelectronic anions. For both the primitive and 
Lennard- Jones models, Equations [14] and [15] show the 
constraints applied to cr, which correspond to the ion 



diameters: 

< cr Li + < (T Na + < cr K + < (T Rb + < (T Cs + 

< o- F - < a cl - < a Bl - < a x - (14) 

^Na+ < °F- 
°K+ < °C1- 
°Rb+ < °Br- 

0Cs+ < (15) 

In accord with the ions' isoelectronic noble gases, the 
Lennard- Jones well depths e also increase going down 
each group (Equation 16): 



< £Li+ < £ Na+ < e K+ < e Rb+ < e Cs+ 

< e F - < e cl - < e Br - < e x - 



(16) 



The scale of the London dispersion interaction, which 
corresponds to the coefficient of the r~~ 6 term in the over- 
all potential, is smaller for cations compared to isoelec- 
tronic anions because of their lower polarizabilities. It 
follows that this coefficient, which is equal to 4eer 6 in 
the Lennard-Jones model, must obey the inequalities in 
Equation [17] 

e Na+( cr Na+) 6 < ^-{^F-f 



e K+(°K+) 6 < eci-^ci-) 6 

eRb+(0Rb+) 6 < eBr-( cr Br-) 6 

£cs+(ct Cs +) 6 < ei-(^i-) 6 



(17) 



Imposition of these constraints focuses the search on the 
subset of parameter space that is physically reasonable 
and promotes transferability of the resulting parameters. 



E. Constrained nonlinear numerical optimization 

We implemented the interaction models and calibra- 
tion objective function with Mathematica 7 (Wolfram 
Research) and used its constrained nonlinear numeri- 
cal optimization routines to simultaneously determine 
all parameters through minimization of the objective 
function. Initial trials showed that for this global op- 
timization problem, the differential evolution methocPS 
achieved better performance than simulated annealing^!, 
the Nelder-Mead simplex methocP^, and local minimiza- 
tion from random initial points. Differential evolution 
is an iterative general-purpose function minimizer that 
evolves a population of solutions to search for the global 
minimum. In our application, each solution is a set of 
parameters for all nine ions. During a single iteration, 
every member of the population competes against a per- 
turbed version of itself for survival. Perturbations con- 
sist of crossing a mutant parameter set with the origi- 
nal such that each parameter randomly inherits its value 
from either the mutant or the original. Mutant sets are 
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generated by randomly selecting three distinct members 
of the population and vectorially adding the first to a 
scaled difference of the other two. If the perturbed so- 
lution improves the objective function score, it replaces 
the original in the population for the next iteration. The 
crossover and mutation processes make differential evo- 
lution robust in the presence of many local minima and 
do not require evaluation of objective function gradients. 

In employing differential evolution, we used a popula- 
tion size of 100, preserved the default options of 0.5 for 
the cross probability and 0.6 for the scaling factor, and 
enabled solution post-processing by the interior point lo- 
cal minimization algorithm^. In both differential evo- 
lution and local minimization, the convergence criterion 
was an absolute or proportional change of less than 10~ 8 
in the appropriate units for all parameters and the ob- 
jective function value. To facilitate convergence of the 
search procedure, we supplied randomly generated initial 
guesses that satisfied most or all of the constraints given 
in Equations [T4ffT7} For each set of initial a guesses, nine 
uniformly and independently distributed random diam- 
eters between 1 A and 5 A were generated, sorted in 
ascending order, and assigned to the nine ions accord- 
ing to a random permutation uniformly selected from 
the fourteen that are consistent with Equations [14] and 



15 For each set of initial Lennard- Jones e guesses, five 
log-uniformly and independently distributed random in- 
teraction energies between 0.001 and 0.75 kcal/mol were 
generated, sorted, and assigned to the five cations, with 
a separate four to the anions, such that the constraints in 
Equation [16] were satisfied but those from Equation |l7| 
were sometimes violated. These bounds on the initial 
guesses did not limit the evolution of the parameters dur- 
ing the minimization process. Initial guesses with con- 
straint violations evolved towards compliance and still 
contributed during the early iterations via mutant gen- 
eration. 

The entire differential evolution process was executed 
100 times with distinct random seeds. Therefore, 10 4 dis- 
tinct initial guesses were evolved in independent groups 
of 100 to produce 100 population champions that we com- 
pared to select the optimal parameter set. In total, each 
model's objective function was evaluated on the order of 
5 x 10 6 times. Despite the robustness of differential evo- 
lution and the large number of parameter sets evaluated, 
the parameter space is high dimensional and global op- 
timality cannot be guaranteed. For the primitive model, 
all 100 populations converged toward one of two solu- 
tions. We selected the one with the better objective 
function score as our final recommended primitive model 
parameter set. In contrast, for the Lennard- Jones model, 
the 100 populations each produced a distinct champion 
parameter set and objective function score, with many 
pushing tightly against the constraints. This suggests 
that the function landscape is rugged with many local 
minima, and that deeper minima corresponding to un- 
physical parameters exist outside the region allowed by 
the constraints. We selected the champion with the best 



objective function score that satisfied all constraints to a 
tolerance of at least 10 -5 in their respective units as our 
final recommended Lennard- Jones parameter set. 



III. RESULTS 

A. Calibrated ion parameters 

Final optimized parameters for both the primitive 
model and Lennard- Jones model are presented in Ta- 
ble [Tl] As expected, all parameters satisfy the periodic 
table trends expressed in Equations [bi}fT7} If satisfac- 
tion of periodic table trends is desired, it is necessary to 
impose the constraints while simultaneously optimizing 
all parameters; calibration protocols that do not include 
these constraints are likely to pr oduce parameters that 
violate them. In some case d 27 * 2 ^ , the attained Lennard- 
Jones e actually reverse the expected trend by decreasing 
down each ion group. In addition to having correct rel- 
ative magnitudes, the ranges of 1.7 to 5.2 A for a and 
0.006 to 0.5 kcal/mol for e are comparable to those of 
Lennard- Jones parameters for noble gasea^, even though 
no absolute bounds other than positivity were enforced 
during their optimization. 



TABLE II. Short-range interaction parameters derived from 
crystal lattice properties. 



Ion 


Primitive model 


Lennard- Jones model 




a (A) 


a (A) 


e (kcal/mol) 


Li+ 


1.716 


1.715 


0.05766 


Na+ 


2.271 


2.497 


0.07826 


K+ 


2.902 


3.184 


0.1183 


Rb+ 


3.165 


3.302 


0.2405 


Cs+ 


3.559 


3.440 


0.5013 


F" 


2.626 


3.954 


0.006465 


cr 


3.600 


4.612 


0.02502 


Br~ 


3.903 


4.812 


0.03596 


r 


4.331 


5.197 


0.04220 



B. Attained values of calibration observables 

The optimized hard sphere diameters for the primi- 
tive model attain a relative root mean square deviation 
(rRMSD) from calibration targets of 4.3%, with RMSDs 
for interionic distances and lattice energies being 0.12 A 
and 8.5 kcal/mol, respectively. While the anions are 
diminished relative to their crystallographically derived 
ionic diameters^, the cations are swollen by a greater 
amount. As shown in Table |TTl} the attained magnitudes 
for all lattice energies and interionic distances are greater 
than their measured values. This indicates that lattice 
energies naively calculated using ionic radii as the hard 
sphere radii would be too negative, necessitating an over- 
all swelling that trades off accuracy in interionic distances 
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in exchange for better accuracy in lattice energies. The 
resulting compromise can be considered a best fit of the 
primitive model to lattice data, and highlights the defi- 
ciency of hard sphere exclusion as a model for short-range 
repulsion between ions. The need to introduce an overall 
swelling of ion diameters to match experimental measure- 
ments has also been encountered in studies of primitive 
model activity coefficients in the mean spherical approx- 
imation^]. 



TABLE III. Negative lattice energies and interionic distances 
attained using lattice-derived primitive model parameters. 
Digits colored blue indicate positive deviations from calibra- 
tion targets. In each number, the most significant colored 
digit is shaded to indicate the magnitude of the deviation, 
with darker shades indicating smaller deviations. 



-LB (kcal/mol) 
ID (A) 


F 


CI 


Br 


I 


Li 


267.3 


218.3 


206.5 


189.5 




2.171 


2.658 


2.810 


3.062 


Na 


237.0 


197.7 


188.0 


175.8 




2.449 


2.936 


3.087 


3.301 


K 


210.0 


178.5 


170.6 


160.5 




2.764 


3.251 


3.402 


3.616 


Rb 


200.4 


171.6 


164.2 


154.8 




2.895 


3.383 


3.534 


3.748 


Cs 


187.7 


163.5 


156.9 


148.4 




3.092 


3.580 


3.731 


3.945 



The optimized Lennard- Jones parameters attain a 
rRMSD from calibration targets of 1.4%, with RMSDs 
for interionic distances and lattice energies being 0.031 A 
and 3.0 kcal/mol, respectively. Table IV shows that these 
deviations are smaller and more balanced than those of 
the primitive model, reflecting an improved ability of the 
Lennard- Jones model to capture the essential physics of 
van der Waals interactions. These deviations are also 
smaller than those obtained by Peng et alP^, which is 
expected since they did not optimize all parameters si- 
multaneously. Compared to Joung and Cheatham's^ 
recommended parameters for use with the TIP3P wa- 
ter model, our deviations are lower for lattice energies 
and higher for interionic distances. However, the sig- 
nificance of these comparisons is limited due to the dis- 
tinct calibration targets used in each study and the vi- 
olation of periodic table trend constraints by the Joung 
and Cheatham parameters. We rejected several candi- 
date parameter sets with better objective function scores 
because they saturated one or more constraints by hav- 



ing multiple ions with nearly identical a, e, and/or r 
coefficients. This suggests that unconstrained minimiza- 
tion of the objective function would lead to even lower 
objective function scores, underscoring the importance of 
the constraints in maintaining physical realism. 



TABLE IV. Negative lattice energies and interionic distances 
attained using lattice-derived Lennard- Jones parameters with 
Lorentz-Berthelot mixing rules. Digits colored blue/red in- 
dicate where attained values exceed/fall short of calibration 
targets. In each number, the most significant colored digit 
is shaded to indicate the magnitude of the deviation, with 
darker shades indicating smaller deviations. 



— IE (kcal/mol) 
ID (A) 


F 


CI 


Br 


I 


Li 


258.9 


210.7 


198.4 


182.5 




2.073 


2.565 


2.732 


2.976 


Na 


226.5 


191.7 


182.9 


170.8 




2.373 


2.822 


2.968 


3.185 


K 


200.7 


172.9 


166.0 


156.6 




2.686 


3.137 


3.278 


3.481 


Rb 


192.7 


166.4 


160.0 


151.2 




2.814 


3.275 


3.417 


3.622 


Cs 


185.4 


158.8 


152.5 


143.8 




2.953 


3.521 


3.675 


3.903 



C. Lattice structure prediction 

Using the optimized Lennard- Jones parameters, we 
compute lattice energies for all twenty salts in both FCC 
and BCC lattice arrangements to check whether the ex- 
perimentally determined crystal structure is correctly fa- 
vored. We find that FCC lattice energies are more neg- 
ative for all twenty salts; this is incorrect for the three 
BCC salts and correct for the other salts. However, the 
gaps between FCC and BCC energies are only 1.6, 1.7, 
and 2.2 kcal/mol for CsCl, CsBr, and Csl, respectively. 
These gaps are narrower than those of all seventeen FCC 
salts, where the correct lattice was favored by between 
2.9 and 22. kcal/mol, and smaller than the deviations of 
attained lattice energies from their experimental values. 
Lingering deviations from calibration targets for the op- 
timized parameters suggests that a limit on the accuracy 
of the Lennard- Jones model has been reached; a more 
realistic model that possibly includes additional param- 
eters is necessary to achieve better agreement with ex- 
periments. This is especially relevant for ions with large 
electron clouds and high polarizability, such as Cs + and 
the larger anions, and may be necessary to correct the 
lattice structure predictions. 



D. Hydration free energies 

To test the transferability of parameters derived us- 
ing lattice properties, we estimate hydration free energies 
AGhyd, an observable that is not used anywhere in their 
derivation. Joung and Chcathanpfil constructed bicubic 
surfaces (a, e) that provide single-ion hydration 

free energies as a function of the ion's a and e by fitting 
the results of several hundred thermodynamic integration 
calculations (note that they used i? m j n = 2( 1 / 6 V instead 
of a). They then combined these surfaces with exper- 
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imental hydration free energies AGS to obtain map- 
pings between a and e such that AG^(tr, e) = AGjjyJ. 
However, the bicubic surfaces themselves, which are not 
influenced by the experimental values, independently 
constitute a useful and valuable tool because they enable 
one to estimate the results of hydration free energy calcu- 
lations for arbitrary Lennard- Jones spheres with mono- 
valent charge to within ~ 0.3 kcal/mol without perform- 
ing any simulations. We employ these surfaces to esti- 
mate the calculated hydration free energy for each alkali 
and halide ion using our Lennard-Jones parameters in 
TIP3PEI, TIP4P-Ew3Sl and SPC/EP water at 298.15 K 
using a standard state of 1 mol/liter for both gas and 
solution phases. As with the original calculations, these 
estimates incorporate Lorentz-Berthelot mixing rules for 
ion-water interactions and omit corrections for the liquid- 
vacuum surface potentiaPHSU. Table [v] compares the 
estimated hydration free energies to experimental mea- 
surements compiled by Schmid et alP^, which also do not 
account for a liquid-vacuum surface potential. 

TABLE V. Negative hydration free energies attained us- 
ing lattice-derived Lennard-Jones parameters and three wa- 
ter models. Experimental measurements from Schmid et al. 
are included for comparison. A temperature of 298.15 K and 
standard state of 1 mol/liter in both the gas and solution 
phase applies to all values. Digits colored blue/red indicate 
where calculated values exceed/fall short of measured values. 
In each number, the most significant colored digit is shaded to 
indicate the magnitude of the deviation, with darker shades 
indicating smaller deviations. 





Calculated —AGhyd 


(kcal/mol) 


Measured — AGhyd 


Ion 


TIP3P 


TIP4P-Ew 


SPC/E 




Li+ 


113.4 


106.1 


112.1 


113.8 


Na+ 


87.6 


82.5 


85.6 


88.7 


K+ 


69.7 


65.8 


67.5 


71.2 


Rb+ 


65.4 


62.1 


63.4 


66.0 


Cs+ 


61.5 


58.9 


59.9 


60.5 


F" 


118.7 


122.4 


123.8 


119.7 


cr 


87.9 


90.0 


90.2 


89.1 


Br~ 


81.5 


83.2 


83.2 


82.7 


I- 


73.4 


74.5 


74.4 


74.3 



The attained RMSDs are 1.0, 4.1, and 2.4 kcal/mol 
for TIP3P, TIP4P-Ew, and SPC/E, respectively. Gradi- 
ents of the bicubic surfaces evaluated at optimized pa- 
rameter values indicate that hydration free energy is es- 
pecially sensitive to Lennard-Jones parameters for small 
ions: across the three water models, the partial deriva- 
tives of — AG^ with respect to a and e are at least 
[31.9,77.2] kcal/mol for Li+ and [38.2,1513.2] kcal/mol 
for F — . Therefore, a deviation of the these parameters 
on the order of 0.1 A for a or 0.01 kcal/mol for e away 
from their lattice-calibrated values could significantly de- 
grade the agreement with experimental hydration free en- 
ergies. The steepness of these bicubic surfaces is intrin- 
sic to the water models and implies that hydration free 
energies are a stringent test of the accuracy of Lennard- 



Jones parameters. Greater discrepancy between simu- 
lation and experiment for TlP4P-Ew compared to the 
three-site models can be explained by displacement of its 
negative charge site relative to the oxygen Lennard-Jones 
center. This increases the distance from the negative 
charge of an ion approaching from the oxygen side and ac- 
counts for the observed trend of increasing severity of hy- 
dration favorability underestimation with decreasing size 
for the cations^. Although caveats regarding assump- 
tions about proton solvation thermodynamics implicit to 
the experimental values apply^, the general agreement 
across water models between simulation and experiment 
on an aqueous system that is completely unlike the crys- 
talline phase used to derive these parameters is evidence 
of their transferability and validity. 



IV. DISCUSSION 

Previous studies have demonstrated the utility of lat- 
tice properties as calibration targets. Peng et alp2l mod- 
ified noble gas parameters by adjusting cation sizes and 
dispersion well depths separately to fit lattice properties 
and periodic table trends within the framework of a 9-6 
van der Waals model with Waldman-Hagler mixing rules. 
Gee et alP*l used lattice dimensions to determine disper- 
sion well depths as part of an effort to reproduce exper- 
imentally derived Kirkwood-Buff integrals using simula- 
tions with SPC/E water. Joung and Cheatham's thor- 
ough study^ fit lattice properties and ion- water binding 
measurements subject to constraints that maintained ac- 
curate hydration free energies, ultimately recommending 
three completely different sets of parameters for three 
water models. They also speculated that it would be 
possible to use lattice properties as the sole calibration 
observables in order to achieve water-independent ion 
parametrization. The present study confirms this specu- 
lation. 

The proliferation of alkali and halide ion parameters for 
nonpolarizable, pairwise additive force fields illustrates 
the inherent disadvantages of solution properties as cali- 
bration targets. Since solvent model parameters are fixed 
during these calibrations, any shortcomings and idiosyn- 
crasies of the solvent model become embedded in the 
resulting parameters. Therefore, every combination of 
solvent model and calibration targets engenders its own 
ion parameters. In fact, since the calculated solution 
properties also depend on simulation design choices such 
as system size, cutoff distance, boundary condition, free 
energy calculation method, electrostatic approximation 
method, etc., careful sensitivity analyses must be per- 
formed during calibration to prevent the resulting pa- 
rameters from becoming specific to those choices as well. 
For small systems, ab initio molecular dynamics and den- 
sity functional theorj^SMI alleviate this problem by ex- 
plicitly modeling electronic degrees of freedom, provid- 
ing a transferable way to simulate ions without deriving 
or choosing a force field. These calculations may serve 



as calibration targets for observables that lack definitive 
experimental measurements. However, the results are 
still sensitive to methodological details such as the choice 
of basis set, pseudopotential, exchange-correlation func- 
tional, etc., and simulations of larger systems such as so- 
lutions of biopolymers at this level of detail are currently 
intractable. 

The ion parameters presented in this study provide an 
alternative to the paradigm of solvent-specific customiza- 
tion for use with nonpolarizable, pairwise additive force 
fields. Since no choices regarding a solvent model or sim- 
ulation design were made during their derivation, they 
are potentially suitable for a broader variety of systems 
compared to ion parameters built on top of a solvent 
model. Lattice-calibrated parameters are immediately 
applicable to salt crystal and molten salt simulations. 
Their independence makes them a reasonable first guess 
in simulations with solvents besides water that lack cus- 
tomized ion parameters. As with water, solvent-ion pair 
potentials determi ned b y the mixing rule may be indi- 
vidually overridde n 1 28 * 58 ! to attain enhanced accurac y for 
solvent-specific properties such as solvation structur e! 59 ! 60 ! 
and gas phase cluster properties^. Retaining the mix- 
ing rule defaults for other pairs preserves the intrinsic 
character of the ions in their interactions with each other 
and with other matter. Encouragingly, the agreement 
with measured hydration free energies shows that true 
transferability, where no such tuning is necessary, is not 
impossible. In fact, these parameters even make it possi- 
ble to reverse the direction of dependence by calibrating 
solvent models while keeping the ion parameters fixed. 
In some situations, aqueous simulations may benefit from 
adopting these parameters despite the abundance of ion 
parameters tuned for particular water models. Lattice- 
calibrated ions are well suited for the ABSINTH implicit 
solvation modeJ221, where ion-water geometries are un- 
available due to the continuum solvent description and 
experimental hydration free energies are direct inputs to 
the model. They are also justified whenever additional 
solutes are present and the balance between ion-water 
and ion-solute interactions is a subject of inquiry; the 
absence of solvent in their derivation makes them inher- 
ently unbiased compared to ion parameters co-derived 
with a water model. 

The general approach of using crystal lattice data to 
derive ion parameters should remain effective for more 
demanding problems such as multivalent ions and polar- 
izable force fields. A wealth of experimental data and 
techniques from solid state physics are already available, 
and the combinatorics that allow formation of many dif- 
ferent salts from few ions readily lead to overdetermined 
systems that are desirable in calibrations. In addition, 
the periodic nature of crystal lattices makes their cal- 
ibration observables inherently easier to calculate com- 
pared to ones involving dense fluids. As demonstrated 
in this study, the use of analytical differentiation and 
polynomial root solvers make calculation of calibration 
observables possible with minimal computational effort 



compared to full equilibrium ensemble simulations, en- 
abling the calculation to be wrapped in a function that is 
called numerous times by a global minimization routine. 
The absence of solvent and symmetry of crystal lattices 
make this calibration process as simple as possible. 



V. CONCLUSION 

Crystal lattice properties are sufficient to determine a 
simultaneous fit of all alkali and halide ion parameters 
for both the primitive model and Lennard- Jones model 
with Lorentz-Berthelot mixing rules. The resulting pa- 
rameters presented here are transferable, consistent with 
periodic table trends, and free from the influence of any 
solvent model. The success and generality of the method 
used to derive them suggests that it may be used as a 
template in future parametrization studies. 
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